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Abstract 

We study the systematic errors of Liischer's formulation of dynamical Wilson quarks and some 
of its variants, in the weak and strong coupling limits, and on a sample of small configurations 
at finite (3. We confirm the existence of an optimal window in the cutoff parameter e, and the 
exponential decrease of the error with the number of boson families. A non-hermitian variant 
improves the approximation further and allows for an odd number of flavors. A simple and 
economical Metropolis test is proposed, which makes the algorithm exact. 

1 Introduction 

Hybrid Monte Carlo (HMC) is now the standard method to simulate fermionic field theories jy]. It 
is exact, relatively simple to implement, and its cost grows just a little faster than the volume V of 
the system [^j. This should not stop us from searching for better algorithms. 

Fermionic interactions give rise, after integration of the fermionic fields, to a non-local determi- 
nant. HMC addresses the problem of the non- locality of this determinant by linearizing the action 
in a succession of molecular dynamics steps. The calculation of the force at each step is obtained 
via the iterative solution of a sparse linear system. But this linearization implies infinitesimal steps, 
which in turn make narrow, very high energy barriers nearly impassable. Thus very long auto- 
correlation times have been observed, eg. for the topological charge Q, and ergodicity may be 
questionable when the vanishing of the fermion determinant divides phase space into disconnected 
regions p]]. A finite step-size algorithm is highly desirable. 

All such algorithms have been abandoned, because of their high cost proportional to V 2 ||. Re- 
cently however, Liischer proposed a finite-step bosonic formulation where one can control, through 
the number of auxiliary bosonic fields, the trade-off between accuracy and computer cost §. Sur- 
prisingly, no systematic investigation of this trade-off has appeared in the literature, except in 
full-blown studies of the complete QCD theory J?], f|. It is our purpose to fill this gap here, by 
studying the weak and strong coupling limits, and the interacting theory on a sample of HMC 
configurations. In Section 2, we recall Liischer's formulation of the QCD determinant, and explain 
our methodology. In Section 3, we present our results for the systematic error of the method. In 
Section 4, we investigate the effect of the usual 'even-odd splitting' of the Dirac operator. In Section 
5, we propose a non-hermitian variant, with improved convergence, suited also for an odd number 
of quark flavors. Section 6 describes a simple Metropolis test which makes the algorithm exact. 
Conclusions follow. 
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2 Liischer's method 



In cases where the fermionic determinant can be written as a square detQ 2 , Liischer has proposed 
a general method to approximate it by a local bosonic action. The essential steps are: 

• find a real polynomial P(x) of even degree n, approaching 1/x as n — ► oo over the spectrum 
of Q 2 , which must be real. 

• decompose P(x) into a product of monomials 

n 

P(x) = c n x n + ... + ClX + c = c n JJ (x - z k ) (1) 

k=l 

such that 

n 

detP(Q 2 ) = c% 11 det(Q - ^)det{Q - JT k ) (2) 
k=i 

where N is the rank of Q 2 . 

• express each factor above as a Gaussian integral over bosonic fields, so that finally 



detQ2 a imw) = ^¥J nwa»*-*> 



hU\rlsh, ] P -<t>i(Q-V^) f (Q-V^k)<t>k 



(3) 



Thus one trades the initial non-local determinant for a bosonic action which is a sum of n local 
terms. This bosonic action can be viewed as a discrete path integral; it converges to the exact 
determinant as the discretization step ~ 1/n tends to zero. Note that the normalization factor c n 
of the polynomial drops out of any expectation value. 

Here we study this approach when applied to dynamical Wilson fermions. Let us denote the 
Dirac matrix by 

D = (1 - kM) (4) 

where M is the lattice Wilson hopping operator and k the hopping parameter. Liischer chooses for 
Q the hermitian operator 

Q = c Q ; Q = l^D (5) 

where cq = (1 + 8k) _1 is introduced to guarantee that the norm of Q is bounded by one, and his 
method to approximate detQ 2 applies to 2 degenerate quark flavors. 

The particular polynomial adopted by Liischer is built from Chebyshev polynomials T n (x), with 
the definitions 

T (x) = 1 (6) 
T 1 (x) = x = cos6 (7) 
T n (x) = cos(n9) (8) 
The approximation polynomial is then defined by 

P(x) = ill^M = Cn f[( x - Zk ) (9) 

x k=i 

where R(x) is the scaled and translated Chebyshev polynomial such that R(0) = 1, 



R(x) = T l 7 l+ l - £j (io) 



with e an adjustable parameter G (0, 1). By straightforward algebra, one verifies that the zeroes of 
P(x), Zk, k = 1, n are given by 
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z k = — — (1-cos — —)-iy/esm — — 11 
2 n + 1 n + 1 

A very important advantage of using Chebyshev polynomials is that the error of the approxi- 
mation can be bounded, and the bound converges exponentially with n: 

\R(x)\ = |1 - xP{x)\ < 2 , Vx G [e, 1] (12) 



The proof is given in section 5. 



Our methodology to assess the error of this method and to compare it with other variants, is 
the following: we compute the quantity 

y = det Q 2 P{Q 2 ) (13) 
by calculating the eigenvalues of Q 2 , Xi,i = 1, ...N, so that 

N 

y = l[X l P(X l ) (14) 
i=l 

To properly cover the full range of gauge couplings, we study the variations of y in 3 different cases. 

i) For free fermions the spectrum of Q is known analytically; we estimate the error by monitoring 

\y 1/N -M- 

ii) In the strong coupling limit we generate a sample of quenched 4 4 configurations, and measure 
the fluctuation of y by monitoring 

1 



a = ^—^v < v 2 >-<y> 2 (15) 

The reason for the normalization factor 1/ < y > is to remove the dependence of the error on the 
arbitrary constant c n (eq.@). 

iii) At intermediate coupling, we perform the same analysis for a set of matrices Q 2 generated by 
hybrid Monte Carlo at j3 = 6 and k = 0.14 on a 4 4 lattice. 



3 Systematic errors of Liischer's method 

In the case of free fermions, our definition {y 1 ^ — 1| of the error depends on the normalization 
coefficient c n . A simple way to compute c n relies on the observation that R{^4^) = 0, so that 

Cn^^fl^-Z,) (16) 
i=l 

We then show in Fig. la the error as a function of the parameter e, for 20, 54, 90 and 148 boson 
fields (k = 0.11, on an 8 4 lattice). As expected, there is an optimal value for e: if e is too large, 
the polynomial approximation degrades over the lower part of the spectrum of Q 2 ; if e is too small, 
the approximation becomes poor over the whole spectrum, because (1 — y / e)/(l approaches 1 

(see eq.(|i~2|)). Not surprisingly, the optimal value of e is slightly larger than the smallest eigenvalue 
of Q 2 , and approaches it as the number of bosonic fields increases. 

At finite (3, fluctuations of \ m in(Q 2 ) complicate the matter: configurations with small X m in{Q ) 
contribute larger errors, so that e must be tuned slightly below < \ m in{Q 2 ) > in order to minimize 



the average error, when this error is small. This behaviour is shown in the (quenched) strong 
coupling in Fig. lb for k = 0.2. It remains unchanged at intermediate couplings (3 = 6 in Fig.lc, 
where we also illustrate a typical spectrum of the matrix D in Fig. 2a. 

From these figures, it appears that the optimal ratio e / X m in(Q 2 ) is relatively insensitive to 0. 
This should simplify the task of tuning e. 

Rescaling the matrix 

As mentioned in |7|], rescaling the matrix Q (see eq.([|)) to Q/cm can be useful. This is equivalent 
to rescaling the roots of the polynomial Zk — > c 2 M z^. In jjj one assumes cm > 1 to guarantee that 
the spectrum is bounded by one. Nonetheless, one may let cm be less than one, if one monitors the 
largest eigenvalue of Q 2 . Under rescaling, this largest eigenvalue approaches 1, making use of the 
full interval [e, 1] of validity of the approximation (|3|). Simultaneously the smallest eigenvalue of Q 2 
increases, with corresponding gains in convergence. One readily sees from equation (|l2^) that the 
number n of boson fields can be multiplied by ~ cm- 

The benefits of tuning cm vanish in the weak coupling limit, since the largest eigenvalue of Q 2 
then tends to 1. Even in the strong coupling, the advantage remains small, O(10 — 20)%. But 
the normalization of Q must be considered carefully when one uses even-odd splitting of the lattice 
sites. 

4 Even-odd splitting 

Liischer's method is based on a polynomial approximation of the inverse, just like iterative methods 
to find the quark propagator, with the difference that the polynomial is predetermined in Liischer's 
case, whereas it is built recursively and adaptively by an iterative linear solver. The similarity of 
the two approaches should encourage us to try here what works well there. In this section we test 
the simplest such idea, that of partitioning the lattice sites into 'even' and 'odd', and of factorizing 
detQ 2 into equal, even and odd factors. 

To justify such a factorization, we simplify the derivation used in || for staggered fermions. 
Define a diagonal operator £ with entries 1 on even sites and —1 on odd sites. This operator 
anticommutes with the hopping matrix M, since M connects even and odd sites. Noting that 
X 2 = 1, one gets 

detD = det(l - kM) = detE(l + kM)T, = det(l + kM) (17) 
Now Z)t = 75D75, and 75 = 1, so that 

detQ 2 = [detD] 2 = det(l - k 2 M 2 ) (18) 

Furthermore we observe that M 2 commutes with E; then if x is an eigenvector of M 2 , so is Y,x, or 
(1 ± S)x, which is non-zero on even or odd sites respectively. Therefore all eigenvalues are two-fold 
degenerate, and one has 

det(l - K 2 M 2 ) even = det(l - k 2 M 2 ) oM (19) 

So we can finally write 

detQ 2 = det(l - K 2 M 2 ) 2 even (20) 

This way we can work with bosonic fields defined on even sites only, saving a factor 2 in memory. 
The bosonic action however is now less local, so that a bosonic update requires about as many 
operations as before. Nonetheless, further gain may come through a reduction in the number of 
bosonic families required for a given accuracy. 
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The reason for this potential gain becomes clear when one looks at the spectrum of a typical 
configuration: the spectrum of Fig. 2a shrinks to that of Fig. 2b after even-odd preconditioning. One 
can now readjust e and vary the number n of bosonic families, with and without even-odd splitting. 
In both cases, the error decreases exponentially with n, as predicted by eq . (|l~2l) ; but it does so much 
faster in the even-odd splitting formulation, as illustrated in Figs. 3a and 3b for (3 = and (3 = 6 
respectively. In both figures the gain is a factor 2 to 3 in n. 

This large gain can be understood by considering the spectral radius of kM, say p = 1 — 7: 
the smallest eigenvalue of entering the determinant eq.(|i~7|) is 7 ^ 1, but that in eq.(|20|) is ~ 27. 
Thus the optimal value of ^Je should also increase by a factor 2, which according to eq.(^) allows 
a reduction of n by the same factor. In addition, one can achieve important gains by tuning the 
rescaling constant cm to bring the largest eigenvalue of Q 2 close to 1, as explained in the previous 
section. The normalization constant which guarantees that the norm of 75(1 — K 2 M 2 ) even is bounded 
by one is now Co = (1 + 64k 2 ) -1 . But with this normalization, the largest eigenvalue of Q 2 becomes 
rather small away from (3 = 00. Thus in Figs. 3a and 3b both, we tuned cm to 0.6 instead of its 
default value of 1. 



5 Non-hermit ian variant 
5.1 Formulation 

Another lesson can be learned from iterative solvers [jlOfl : the method of biconjugate gradients 
(BiCG) [ll]] requires less work than that of conjugate gradients (CG) fl^ ]. BiCG applies Lanczos 
polynomials on the matrices D and to construct the solution, whereas CG applies the Lanczos 
polynomial on the matrix Q 2 . This observation motivated us to look for an approximation to 
1/D itself, instead of l/Q 2 . This should be possible as long as detD is positive. But a negative 
determinant can only be caused by negative real eigenvalues; this situation can only occur for very 



small quark masses as studied in the quenched case by [13|, beyond those reachable by present-day 
simulations of full QCD. The benefits of the approach proposed below are two-fold: convergence of 
the approximation is improved; and the simulation of an odd number of quark flavors (or of any 
number of non-degenerate flavors) becomes possible. 

Let P{z) be a polynomial of even degree n, with real coefficients, defined in the complex plane 
(since the spectrum of D is complex), with complex conjugate roots z k ,Imz k ^ 0, k = l,...,n. 
Assume that this polynomial satisfies P(z) — > 1/z for n —* 00, Vz G S, where S is a domain in the 
complex plane containing the spectrum of D. Then one can write 

n/2 

detP(D) = c% det(D - z k )det(D - z k ) (21) 

k=l 

Using the fact that D = 75.0^75 and 75 = 1 one gets 

det(D - z k ) = det{D ] - z k ) (22) 



It follows that 



n/2 

detP(D) = c n Y[ det(D - z k )Uet(D - z k ) (23) 

k=l 



and, in analogy with eq.(||) 



n/2 



k=l 



The approximation polynomial can be defined as before by 



P(z) = ^ = c n J] (z - z k ) (25) 

z k=l 

where R(z) is a polynomial of degree n+1 such that R(0) = 1. |i2(z)| should be small in the domain 
S of the approximation. The spectrum of D has a more or less elliptical shape, as can be seen from 
Fig. 2a. It is then natural to choose for the boundary dS an ellipse. In that case (provided the origin 
is not inside the ellipse), it turns out that the polynomial which minimizes maxs\R{z)\ is again a 
Chebyshev polynomial ||l4j| . 

The same definitions eqs.(|6||8|) remain valid for T n (w), where now w = x + iy and 9 = 9\ + 162, 
that is 



w = cosO = cosOi COSI162 — i sin6\ sinh02 (26) 

where 6\ G [— tt, tt] and, for instance, 62 € [0, +00 [. Note that, for 62 fixed, w(9%) describes an ellipse 
in the complex plane, with foci ±1 and semiaxes cosh02 and sinh02- 

From these definitions we can now specify the form of R(z). Choose for dS the ellipse centered 
at (d, 0), with large semiaxis a < d, and focal distance c < a. Then the scaled and translated 
Chebyshev polynomial R(z) is defined by 

rp I Z — d \ 

m = T n+1 c d (27) 

From this definition one can show, in analogy with the hermitian case, that the roots of P(z), 
Zk, k = l, n are 

2;^ = d(l — cos —) — i\J d 2 — c 2 sin — (28) 

n+1 n+1 

so that the z^'s lie on the ellipse of same center and foci as dS, which goes through the origin. 

We can now derive an error bound for the above polynomial approximation in the complex 
plane. It states that 

1^)1 = H-™ 2 .V, £5 (29) 



Proof. By definition 



where 



= cos(n + l)g 

v ' cos(n + 1)9 v ' 

z — d d 
cos 9 = , cos#o = — (31) 

c c 



From eq. (|26|) , one sees that 9q = it + ia, where a = cosh l d/c. Writing 9 = 9\ + i92, one gets 



R( z ) = — — — —(cos(n + l)6»i cosh(n + 1)0 2 - % sin(n + 1)9 X sinh(n + 1)9 2 ) (32) 

cosh(n + lja 

so that 

< cosh(n+l)e 2 
1 Wl ~ cosh{n + l)a 

The coordinates of any point z in S can be expressed as (^i,^). Successive 6Vs define nested 
ellipses, all of center d and focal distance c: the innermost one, for 02 = 0, is the real segment 
[d — c, d + c]; the outermost one, for 62 = 9 ma x = cosh -1 a /c, is dS. Along each such ellipse |-R(^)| 



is bounded by C cosh{n+i)a > an< ^ reacnes this bound when z is real. The bound increases with 62, so 
that over S, \R(z)\ is bounded by 

cosh(n + l)9 max = ein+memax -a) 1 + e^+^ 
cosh(n + l)a 1 + e -2(n+l)a ^ ^ 

The second factor is bounded by 2. Substituting 9 m ax and a by their definitions in terms of 
a,c,d, one recovers eq. (|29|) . 

This derivation makes it clear that the error |-R(.z)| is maximum on the real points of dS, but 
decreases exponentially inside S. This is in sharp contrast with the hermitian approximation of 
Liischer. 

5.2 Optimization and comparison with the hermitian case 

Two limiting cases for dS are of special interest: the circle (c = 0), which corresponds to the 
spectrum of D in the (quenched) strong coupling; and the real line segment [d — a,d + a] (c = a), 
which corresponds to the spectrum of Q 2 originally considered by Liischer. For these 2 cases the 
error bound eq.(|2^) becomes 2(|) n+1 and 2( d+ ^° 2 a z ) n+l respectively. To make contact with the 
hermitian error bound eq.(^), we just express a = and d = We recover then eq.(^) over 
the real segment [e, 1]; over the circle we get: 



1 



n+l 



1^)1 ^ 2 { — ) (35) 

We are now in a position to compare the hermitian and non-hermitian approximations, in the strong 
coupling limit. Call 77 ~ m quar k the smallest eigenvalue of D. Then the smallest eigenvalue of Q 2 
will be rj 2 , and the hermitian and non-hermitian bounds eqs.(|T2|) and ( |35| ) are identical^]. It is true 
that only n/2 bosonic fields are necessary to generate det(D) in the non-hermitian case; but for 2 
degenerate quark flavors, one needs another n/2 fields to finally approximate det(D) 2 , for the same 
total of n fields as in the hermitian case. 

When (3 increases, the error bound for the non-hermitian case improves. The spectrum of D can 
then be contained in a more elongated ellipse whose aspect ratio tends to 2 in the free field limit. 
It is straightforward to verify from eq.(|29|) that the convergence rate is multiplied by 2, allowing a 
reduction of n by the same factor. 

As importantly, the error in the non-hermitian approximation decreases exponentially inside S. 
To exhibit this difference with the hermitian approximation, we show in Fig. 4 the magnitude of the 
error along the real axis, for r] = 0.1 and n = 20. That is, the range of the approximation in the 
hermitian case is [0.01, 1], whereas in the non-hermitian case it is [0.1, 1]. In the hermitian case, the 
error oscillates with constant amplitude 0(1O~ 2 ) over the interval [rj 2 , 1]. In the non-hermitian case, 
when dS is a circle the error at x = r\ and x = 1 is also 0(1O~ 2 ); when dS is an ellipse of aspect 
ratio 2, the error at x = r\ and x = 1 is about squared. Either way, the error falls off exponentially 
inside the approximation range, until it becomes uniformly oscillating in the segment [d — c, d + c] . 
So the accuracy on interior eigenvalues of D is much improved. 

We make numerical tests of this non-hermitian approximation for HMC and strong coupling 
configurations. We apply the same methodology as before, but in the non-hermitian case we measure 
the fluctuations of the quantity 

N 

y = l[(X t P(X l )) 2 (36) 



1 We greatly simplify here a subtle issue: if the smallest singular value of D is rj, then the smallest eigenvalue of 
Q 2 is rj 2 ; but the relationship between eigenvalues of D and of Q 2 could be quite different, especially in the chirally 
broken phase. In particular one easily proves that a m in{D) < \X\min{D), so the comparison we make here is a worst 
case scenario: the non-hermitian variant will always be at least as accurate as the hermitian one. 
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where now Aj,i = 1,...,N are eigenvalues of the quark matrix D. With this definition of y we 
can compare the results directly to the hermitian case. In Figs. 5 we show that the errors of the 
non-hermitian approximation are reduced faster than in the hermitian approximation. From Fig. 5a 
we can see that in the quenched strong coupling, the gain is a factor 4 to 5, increasing as the quark 
mass is reduced. This large gain is caused by the high density of interior eigenvalues. In Fig. 5b 
at (3 = 6, the gain is a factor ~ 1.5, coming mostly from the elliptical shape of dS. Note that 
simulating one flavor to the same accuracy would require half as many fields. 



6 Metropolis test 

Liischer's original proposal includes the monitoring of the error, and the possibility of obtaining 
exact results by re-weighting the Monte Carlo measurements of each observable. The natural way 
to calculate the error eq.([D|) is to express it in an eigenbasis of Q 2 , obtained by the Lanczos 
algorithm. However two obstacles appear: the Lanczos algorithm is affected by roundoff errors, 
which can be controlled only if eigenvalue multiplicities are known 111], and its cost grows like the 
square of the volume V of the lattice, so that it becomes overwhelmingly expensive on large lattices. 
Nevertheless attempts at using the Lanczos method in a Metropolis test show that most if not all 
of the error can be removed [|T6| , [|. 

In fact, it is sufficient to construct an unbiased estimator of the ratio of errors between the new 



and the old configurations. One can then use the noisy Monte Carlo method of [17], which was 
successfully applied to fermionic simulations before the advent of Hybrid Monte Carlo [||, 18 1. 

The ratio to estimate is det(D' P(D')) 2 /det{DP(D)) 2 , calling D' and D the Dirac operators for 
the new and old configurations respectively. Taking the denominator as a partition function, this 
ratio can be rewritten 

< e -n^W-i)v > (37 ) 

where the average <> is taken over all Gaussian vectors r/, and W = [D' P(D')}^ 1 DP(D). It is 
sufficient to estimate this ratio by taking one Gaussian rj only. This requires the solution of a linear 
system, which will take just a few iterations since the matrix D'P(D') is almost the identity. The 
cost of this additional step is therefore similar to that of an update of the 4>k's. As the volume V of 
the lattice grows, one should keep constant the acceptance of this Metropolis test. This is achieved 
if the error per eigenvalue scales like V . Because the approximation is exponential in n, this can 
be accomplished by a modest increase oc LogV in n, and a CPU cost scaling like V(LogV) 2 . Under 
this rescaling of V, the number of iterations necessary to solve the linear system above remains 
constant, so that the overhead of the Metropolis test, measured in update sweeps, does not change. 
Similarly, when the quark mass is decreased but the acceptance is kept constant, the overhead of 
the Metropolis step, measured in update sweeps, will not change. 



7 Conclusion 

We have studied the systematic errors of Liischer's method to simulate dynamical quarks, in the full 
range of strong to weak coupling. Two parameters define the approximation, e and n: e is the lower 
cutoff of the polynomial approximation, n is the number of auxiliary bosonic families. We confirm 
the importance of tuning e near the minimum eigenvalue of Q 2 ; the optimal ratio e/ < \ m in(Q 2 ) > 
approaches 1 from above as n increases, in a manner almost independent of /?; it actually becomes 
slightly less than 1 for large n, because configurations with small \ m in{Q 2 ) give larger errors. For 
e fixed, we confirm the exponential convergence of the approximation with n. 

In addition, we have studied the improvement of even-odd preconditioning: it allows a reduction 
of n by a factor 2 to 3. A full reduction however can only be achieved by tuning to less than 1 the 
normalization parameter cm of the matrix Q. 

We have introduced a non-hermitian variant, which allows the simulation of an odd number 
of quark flavors. The number of bosonic families is proportional to that of flavors. For 2 flavors, 



the non-hermitian formulation allows a reduction of n by a factor 1.5 to 5. Full Monte Carlo tests 
of this variant are in progress. They will shed some light on the critical dynamics of the bosonic 
fields 4>k, which depend on the singular values of {D — z^), and might be different from the original 
hermitian version. 

Even-odd preconditioning and non-hermitian formulation are most advantageous for large quark 
masses and small (3 respectively. These two improvements can be combined when appropriate. 

Finally we have proposed a Metropolis test which removes any approximation. The overhead of 
this Metropolis test is modest and independent of the lattice volume and the quark mass. 

The error we have considered is that on the fermionic determinant itself. The error measured 
on a given observable during a Monte Carlo simulation of the Liischer action will depend on the 
overlap of that observable with the various eigenmodes of the determinant, and may be smaller than 
we measured here. 
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Figures: 



• Figure 1: the magnitude of the error on the determinant as a function of e, for n = 20, 54, 90, 148, 
in the case of free field (Fig. la), (quenched) strong coupling (Fig. lb), and at (3 = 6 (Fig.lc). 
k = 0.11, 0.2, 0.14 respectively. Arrows mark the minimum and the average of X m in(Q 2 ) over 
the ensemble of configurations. 

• Figure 2: a typical spectrum of the Dirac matrix D, eq.(||), in the complex plane (4 4 lattice, 
P = 6, k = 0.14). In Fig. 2b the spectrum is shown after even-odd preconditioning. 

• Figure 3: the magnitude of the error on the determinant as a function of n, in the original 
formulation (dotted line) and after even-odd preconditioning (solid line). Fig. 3a corresponds 
to = 0, k = 0.2, Fig.3b to P = 6, k = 0.14. 

• Figure 4: the magnitude of the error of the polynomial approximation P(x) as x varies 
from to 1. The three approximations shown are the original hermitian approximation of 
Liischer, the non- hermitian approximation inside a circle as appropriate for (3 = 0, and the 
non-hermitian approximation inside an ellipse of aspect ratio 2, as appropriate for (3 = oo. 
All cases correspond to the same quark mass. 

• Figure 5: the magnitude of the error on the determinant as a function of n, in the original 
formulation (dotted line) and in the non-hermitian variant (solid line). Fig. 5a corresponds to 
P = 0, Fig.5b top = 6. 
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